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ABSTRACT 

A numerical method for two-dimensional cavitating flow 
is developed for the flat plate. The linearized boundary 
value problem is restated as a set of coupled integral equation*. 
The integral equations are approximated numerically. The 
numerical approximation is executed by a Fortran IV computer 
program. The computed results are compared to the analytic 
solution. This method should provide insight into developing a 
method for three-dimensional cavitating flows and is readily 
extendable to cambered profiles. 



Thesis Supervisor! Patrick Leehey 

Titlei Professor of Applied Mechanics 
Professor of Naval Architecture 



ACKNOWLEDGEMENTS 



The guidance and invaluable constructive criticism 
provided by Professor Patrick Leehey is gratefully acknowl* 
edged. The very useful suggestions of Professor Justin 
Kerwin are also gratefully acknowledged. 



TABLE OP CONTENTS 



4 



Page 

ABSTRACT 

ACKNOWLEDGEMENTS 3 

LIST OP FIGURES 5 

LIST OF TABLES 6 

NOMENCLATURE 7 

CHAPTER I INTRODUCTION 9 

CHAPTER II THE LINEARIZED BOUNDARY VALUE PROBLEM 12 

CHAPTER III THE INTEGRAL EQUATION FORMULATION 18 

CHAPTER IV NUMERICAL APPROXIMATION TO 

THE INTEGRAL EQUATIONS 24 

CHAPTER V ANALYTIC SOLUTION 34 

CHAPTER VI RESULTS, CONCLUSIONS 

AND RECOMMENDATIONS 39 

REFERENCES 66 

APPENDIX A COMPUTER PROGRAM 68 

APPENDIX B SAMPLE COMPUTER OUTPUT 88 

APPENDIX C COMPUTER PROGRAM FOR COMPUTATION OF 

CAVITY LENGTH FROM ARBITRARY VALUES OF 

ANGLE OF ATTACK AND CAVITATION NUMBER 94 



LIST OP FIGURES 



5 



Figure No. 




Page 


1. 


Coordinate Systems 


45 


2. 


Element Arrangements 


46 


3. 


Variation of Last Source 


47 


4. 


Mapped Planes 


48 


5. 


Vortex Distribution, Jl * o.5 


49 


6. 


Source Distribution, 0.5 


50 


7. 


Cavity Length - Partial Cavitation 


51 


8. 


Lift Coefficient - Partial Cavitation 


52 


9. 


Cavity Area » Partial Cavitation 


53 


10. 


Vortex Distribution, Jl * 1.5 


55 


11. 


Source Distribution, Jl * 1.5 


56 


12. 


Cavity Length - Super Cavitation 


57 


13. 


Lift Coefficient - Super Cavitation 


58 


14. 


Cavity Area - Super Cavitation 


59 


15. 


Numerical Convergence } sl-o.s 


61 


16. 


Computation Time 


62 



LIST OP TABLES 

Page 

Values of M and N for 54 

figures 5 through 9 

Values of M and N for 

figures 10 through 14 60 

Computed Results 63 

Singular Behavior 65 

Results of Cavity 

Length Computations 97 



NOMENCLATURE 



A = Cavity Area 

c = chord length 

[Cl = coefficient matrix 

c j L = element of coefficient matrix 

Cp = pressure coefficient 

Cl a lift coefficient 

f = location of vortex control points 

f - location of source control points 

h = locus of points describing mean camber line 
and cavity surface 
Z * cavity length 

M * number of vortex elements 

N = number of source elements 

P = pressure 

P*, “ pressure at infinity 

q = velocity 

q(.p = source distribution 
u,v = perturbation velocity components 
= inflow velocity 
w = complex velocity = u-iv 
x,y = coordinate system 

x^ * lower boundary of vortex or source element 
Z 



complex plane * x + iy 
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(“*) = quantity with a bar is a dimensional quantity 

( ) = no bar indicates a nondimensional variable 

( ) + * variable to be evaluated at x = 0 + 

( )“ * variable to be evaluated at x * 0” 

( = variable to be differentiated wrt x 

(X * angle of attack 

= vortex distribution 
£ * mapped plane + 




y[ = dummy coordinate system for integration 




P « fluid density 



(5 * cavitation number 
= velocity potential 
CD = perturbation velocity potential 



n = degree of singularity 
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CHAPTER I 
INTRODUCTION 

The objective of this investigation is to develop a 
numerical solution to cavitating flows which can be compared to 
a known analytic solution. Therefore, the cavitating two-dimen- 
sional flat plate in steady flow has been chosen for this 
investigation. Geurst [l,2,3j has given solutions for the 
linearized problem for both partial and super cavitation. 

The purpose in considering this problem is to determine 
and overcome the numerical difficulties associated with this 
problem. Adequately solving the two-dimensional problem is a 
preliminary step in developing a numerical solution for the far 
more difficult three-dimensional and unsteady flow problems. 

Purther the numerical method can be applied to foils 
of different camber lines. Geurst [2] has obtained a camber 
line solution for the case of parabolic camber in partial 
cavitation only. 

Widnall [k] developed a numerical method for three- 
dimensional, unsteady, supercavitating flows. However, Widnall 
used assumed cavity lengths and demonstrated that, for long 
cavities, lift coefficients were relatively insensitive to 
cavity length. Clearly for short super cavities, on the order 
of chord length, and partial cavities, a proper statement of 
cavity termination is required for the proper solution of 
cavitating flows. 
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In this respect this report is a preliminary investiga- 
tion into the inclusion of cavity termination in cavitating 
flows. The results of this report should be applicable to 
three-dimensional and unsteady cavitating flows. 

The solution of two-dimensional, cavitating flows pre- 
sented herein envolves five distinct steps. These are i 

1 ) statements of the linearized boundary 
value problem, 

2) solution of the boundary value problem in 
the form of coupled integral equations in 
source and vortex distributions, 

3) a numerical approximation to the integral 
equations, 

4) execution of the numerical approximation via 
FORTRAN computer program, 

and 

5) validation of the numerical results. 
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Linearizing the boundary value problem and rewriting 
the problem in the form of integral equations is not the only 
method of solution. Another possible approach is the finite 
difference method. Jeppson [5»6j and Mogel and Street |V] have 
taken the finite difference method as far as a circular disk 
perpendicular to the main stream with an infinite cavity. The 
major disadvantage of the finite difference method is the large 
programing and computational effort required. 

The integral equation method for the linearized problem 
represents a far smaller computational effort. It is for this 
reason, plus the existence of the analytic solution, that the 
integral equation method is used. For an excellent discussion 
of the various methods of solution see Birkhoff [8] • 

The following chapters of this report will discuss the 



stated steps of the solution. 



CHAPTER II 
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THE LINEARIZED BOUNDARY VALUE PROBLEM 



Geurst £1,2,3} has presented a development of the 
linearized boundary value problem for both partial and super 
cavitation. A similar development is given here, except that 
the coordinate system and the nondimensionalization of variables 
differ. These differences are strictly a matter of convenience. 
Geurst *s method of solution relies on a series expansion of the 
complex velocity in the neighborhood of infinity. This in 
turn requires that the boundary value relations to local slope 
on the foil and cavitation number be stated at infinity. 

The neighborhood of infinity is not convenient to 
computer solutions. Therefore the boundary value relations to 
local slope on the foil and cavitation number need to be stated 
on the foil and cavity. 

The coordinate system and the foil-cavity relation to 
the coordinate system are shown in figure 1(a). The inflow 
velocity, , is assumed to be parallel to the x axis and 

uniform in the y and z directions. The angle of attack, o(. » 

is taken to be the nose-tail line for cambered profiles. 

First the steady flow form of Bernoulli’s equation is 
written between infinity and another point in the flow* 
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This equation is now rewritten in the form of a nondimensional 
pressure coefficient. 

<» =(*>€<: -i). 

A \ °o 

On the cavity the pressure coefficient is the cavitation 
number, (J ~ • The velocity potential can be written as the sum of 
the potential due to the inflow velocity plus a perturbation 
potential* 

O) &Z + 0 



Where 



(£) * perturbation velocity potential < 



Prom equation (3) the x and y components of velocity are* 



(4) 



*L-0 + if = o + a 
lx ~ 3* 



C5) 



ii- l£- v 

>9 >9 
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Where 

u * x component of the perturbation velocity 

v = y component of the perturbation velocity. 

It is assumed that the perturbation velocity is very much less 
than the inflow velocity, that isi 

U, « Voo 

On the cavity we now have* 

+- x “cQo *££ + 

using(6) 

w M = / + z *c/a„ 

Thus equation (2) for the cavity becomes t 

(8) u. c = % 

Where 

u c= a ‘/Q„ 



and Uco is now used as a nondimensionalization constant. 
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Equation (8) is the statement of the linearised boundary con- 
dition on the cavity surface* The next step is to obtain the 
boundary condition on the wetted surface of the foil* 

On the foil wetted surface the boundary condition is 
flow tangency. 

-y~ _ cJ A 

(9) 0„ +a ~ d* 



Which becomes 



or 

CIO) 



— 7 — -h +- higher order iersns - —— 

LU * <** 







Also, it is assumed that the foil is very thin compared 
to the chord length such that the boundary condition in equa- 
tion (10) can be applied to the mean camber line. 

Further the entire foil cavity system is collapsed to 
the x-z plane and the boundary conditions applied on the x-z 
plane with all length dimensions nondimensionalized by the chord 
length, figure 1(b)* Another condition to be satisfied is cavity 
termination. 

Geurst ^2j proved that for the linearized problem, the 
re-entrant jet and Riabouchinsky models for cavity termination 
reduce to a statement that the cavity closes at its end. This 
condition can be written as 



The next condition to be satisfied is the condition at 
the trailing edge of the foil (x * 1). This condition requires 
that there be no pressure jump accross the foil at the trailing 
edge. 

(is Cp (i) - Cp(o =o 



In summary the conditions that the solution must 



satisfy 


are t 




(c/) 


u* = <r/z 


0 4 * -J? 


(Cl) 


v* h x 


on wetted surface 
of the foil 


(C3) 


\\ k; - d = o 




(Cl) 


c.pO) -c;o)=o 
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The last condition is that the perturbation potential 
satisfies Laplace* s equation i 

(C5-) v a j2 9-o. 

The next chapter will give a solution to Laplace's 
equation which satisfies the above boundary conditions* 
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CHAPTER III 



THE INTEGRAL EQUATION FORMULATION 

This chapter presents the development of the integral 
equation formulation of the problem* This development proceeds 
from the solution for the velocity field induced by a two- 
dimensional distribution of vorticies and sources in space. 

Then it can be shown that this velocity field will satisfy the 
boundary conditions. This results in a set of coupled integral 
equations. 



First consider the flow to be in the complex @ plane. 



Where i 

2 ~ X + ly 

and the complex velocity is given byi 

urcz) - u — i'ir 

Then the appropriate general form of the solution is t 





or 
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(15) w 



«) = skT 

Jo 



where , 



t 

i 






s t> 



<g(f)= source distribution over the cavity 
YCg) m vortex distribution over the foil 

1 for a partial cavity 
a =] 

for a super cavity 

0 for for a partial cavity 



and 



U$) * 0 for ± < J 



for a super cavity. 



Plemelj's formulas can be applied to equation (15) 
as follows £l 03 * 



Cl 6) \ai + ( x) — w~ (x) = - YCx) — i%C^) 

or) w'm * jg 
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When equations (16) and (17) are added together and 
the real part is taken, the result isi 

rX 

«'<*> - - * ih f 

"O 

or, with the cavity boundary condition, 

(IS) <r - - Y(x) -i - 7?j Jg 

with 



XW'O for PC 2./ 



The integral in equation (18) is a Cauchy principal value in- 
tegral over the length of the cavity* 

When equations (16) and (17) are subtracted from each 
other and the imaginary part is taken the result ist 



V-'cx) - - 1C5-? + -L 

v (X) - A + 2^ 



mi 

*-s 






or, with the wetted surface boundary condition, 

us) hi ex) = - + akf j s 



with 



<%CX) = O for X >J 

Again the integral in equation (19) is a Cauchy principal 
value integral* For a flat plate at an angle of attack 

(3.0) fj x) = - c* 



21 



and this is the case taken throughout the rest of this report. 
The source distribution q(x) represents the slope difference 
between the cavity surface and the mean camber line of the 
foil at x. While the vortex distribution Y (x) represents the 
difference in the x component of perturbation velocity between 
the upper and lower surface at x. 

(AO <gpO = i r+(x) - ns~(z) 

(iU) y (m) - -u + (x) + u~ (Z) 

These statements follow directly from equations (16) and ( 17) • 
When equations (2), (7) and (22) are combined the jump in 
pressure on the foil isi 



The combination of condition (C 3) and equation (27) gives, for 
the closure condition, 



It is clear from equation (23) that condition (C 4) is* 



(Xi) C* w -Cf(-x) = -ZYC -X) 



and the coefficient is* 





X(l) = o 
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Equations (18), (19) » (20), (25) and (26) constitute a complete 
problem statement for the steady flow cavitation of a two- 
dimensional flat plate, only, it is convenient to rearrange 
equations (18) and (19) into a form that contains the ratio 
of angle of attack to cavitation number. 




Equations (27), (28) and (29) » subject to Y(J) = o , are 
sufficient to obtain cavity length, the source distribution 
and the vortex distribution as functions of the ratio ot/<r only. 
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The next section will develop the numerical approxima- 
tion to equations (2?) t (28) and (29) and the method of solution* 
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CHAPTER IV 

NUMERICAL APPROXIMATION TO THE INTEGRAL EQUATIONS 

The method used in the numerical approximation is to 
approximate the integral equations (27)* (28) and (29) by a 
finite set of linear algebraic equations with N unknown source 
densities, M unknown vortex densities and the ratio oc/o- • 

The numerical integration procedure requires that cavity length 
be the independent variable and be computed. 

This method requires an assumption about the functional 
form of the source and vortex distributions. It is assumed that 
the vortex and source distributions are constant over a small 
element of the foil or cavity, figure 2. The last vortex 
element is assumed to be linear with a value of sero at the 
trailing edge to satisfy equation (26). The amplitude of each 
1 1 is unknown. 

In each vortex element a control point is placed at 
which the flow tangency boundary condition is satisfied. In 
each source element a control point is placed at which the 
cavity boundary condition is satisfied. The vortex control points 
are at iXj and the source control points are at • The 

relations for the control point locations are* 

1 ) for the vortex control 



(30) 



-- *4 + (X,J„ - %j)f 5 



e>< / 



25 



and 



2) for the sources 



(30 -*,j - Xjij + ( ~Xjj)f s °<f s < I 



The fractions f and f 5 are yet to be determined on the 
basis of obtaining stable results. 

There will be one equation for each vortex control 
point, giving M equations. There will be one equation for 
each source control point, giving N equations. Plus there will 
be one equation for the cavity closure condition. The vortex 
equations are t 




and for the M^ h control point 




+ 




Where , 



■= O for j > A' 



hr 






_ 



The source equations are t 
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(34) 71 — 71 (j- 



r%ti+ 1 

*■ + 



Xjj+I 



+' 



'*1 



I 

. *j - f 3 %-° for j > 



M 



and the closure condition 



A/ 



■ Xjii+i 



( 35 ) 



fer J *£ 



_iL 

0" 



<4* = 6? 



Equations (32) through (35) are the numerical approximation to 
the integral equations* When the indicated integration is per- 
formed equations ( 32 ), (33) » (3*0 and (35) become 



(36) 



o=.a.v$ - «# + %rM-ih) + 



i=l 



+ ¥l~hlL JU ’(^ L ) +0 , #*of.r j>" 



(37) o =ang -ir p + %*[/ + 



< = / 



Xai \ 

Xa\ Jiei+t ) 



0 




o ■ (or 



+ 

Al > N 



2 ? 



(38) if - -V 



* * **■<&> +£***#*) 






and 



«• _ 



cr 



— O for j > N 



m 



N 







“ (x M „ - X M ) = O 



L=! 



Equations ( 36 ) through (39) are linear algebraic equations in 
the unknowns -§p , and -~r • It is convenient to rewrits 
these equations in matrix form. 
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The coefficient matrix, [C], is formed from the coeffi- 
cients of £ , li f JL , M- and in equations ( 36 ) 
through (39). 

The first submatrix, (I), is formed as follows by the 
coefficients of r in equations ( 36 ) and (3?)* 

£ 

W) c j j - ^ , J * m 

r^) c„„ = /+ 

(43) Cji - £y> ( Tj _ ) 0 c t A7 j ? /v, 

W ' /- ^ ^ f 4- / J j* A/ 

The second submatrix, (II), is given by the coefficient 
of ^/c" ®<1 unions (36) and (37)i 

( 4 5) = “ ^ 

(46) C ji - O 5 M+! £ L ± Mm 

an d / 9^ j + Ai . 
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Submatrix (III) is the coefficient of <x/cr in equations 



06) and (37) « 

(4-7) Cj )AA +H H = 27T for I±j£M 

(40) Cj )M +hh =0 for A1 <-j± M+A/^/ . 



The source equations are contained submatricies (IV), 
(V), and (VI), The submatrix (IV) contains the coefficients 
of %rin equation (38) i 

<T49J = -7T, 



and 



(50) C H = o 



other wi s& ■ 
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Submatrix (V) contains the coefficients of JL in 

<r 

equation (38)* 



(5/) 



C 






= JL~( 



JH 



xj 



- ^t>/ / } 



A 1 +• j ^ /*f +£ 
/ ^ jj i 6 A/ 



and 

(5<2) 






The last submatrix (VI), contains the closure condition, 
equation (39)* 



( 53 ) (-MtM+iji. 



X/Uti ~ *X-jU i 

o 



M + > ± L ± M + W 
otherwise 



To solve equation (40) for the unknown vortex and source 
distributions and the ratio <*Ar the coefficient matrix is in- 
verted. The inverted coefficient matrix is then multiplied by 
the boundary condition matrix. 

Thus the solution can be written as* 



(S+J 







<r 




fH 




<r 
%\k r 


— 


H 








ZT 






o 



o 

It 




On the surface this procedure appears to be straight- 
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forward and should produce a solution without difficulty. 
However, this is not the case. 

Difficulties in obtaining a meaningful numerical solution 
lie in the selection of control point locations within each 
source or vortex element. Control point locations are critical 
to approximating the singular behavior of the analytic solution. 
It was found that vortex control points should be at 90 percent 
of the element length and source control points at the element 
midpoint. However, special care must be taken with a few of the 
source control points. 

The critical control points are the first and last 
source control points. A systematic variation of the control 
point position determined that the last source control point 
should be placed at 10 percent of the element length. As an 
example of the importance of this control point, the variation 
of the computed ratio of with control point location is 
shown in figure 3 for a cavity 0.5 times chord length. 

The selected control point position is at 10 percent 
of the element length. Variation of the first source control 
point indicated the placement should be at 90 percent of the 
element length. Since the source distribution is singular at 
the cavity leading edge and termination, these control points 
should be placed away from the singularities. 



Variation of other control points (second and next to 
last source, last vortex, first source past trailing edge for 
super cavity) showed the best placement to be the generalized 
locations* Thus the fractions f and f 8 arei 

f ®0*90 for all vortex elements, 

f s * 0.90 for first source control point, 

f s * 0.10 for last source control point, 

and 

f s * 0.50 for all other source control points. 

The computer program that executes the solution to the 
described numerical procedure is listed in Appendix A. This 
program is composed of a main program, which performs Input/ 
Output operations and logic control, and three subroutines. 

The first subroutine, CPGEN, determines the vortex and 
source element sizes and the locations of the control points. 

The required input information for this subroutine is the cavity 
length, the number of vortex elements and the number of source 
elements. 

The second subroutine, MATRIX, computes the coefficient 
matrix, iC ], and the boundary condition vector. The required 
input information for MATRIX is the output of CPGEN and the 
number of vortex and source elements. 

The last subroutine, RMINV, inverts the coefficient 
matrix. This subroutine is a standard MIT matrix inversion 
routine [ll] • The inverted matrix is then multiplied by the 
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boundary condition vector, in the main program, to obtain 
the solution vector for and -J? . 

A special note on the program notation is made here* 
The vector XU(I) is the upper boundary of a vortex or source 
element. This vector is given byi 

XU(I) =? XL( I + 1) 

through an assignment statement in CPGEN* 



CHAPTER V 
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ANALYTIC SOLUTION 



The method of gauging the numerical results is the 
analytic solution for the cavitating flow of two-dimensional 
flat plate. The first half of this chapter will discuss Geurst*s 
Solution for partial cavitation £1,2] . The second half will 
discuss Geurst's solution for super cavitation £ 3 J • 

Geurst solves the partial cavitation problem by a 
conformal mapping of the *2 plane to the £ plane as follows. 
First the % plane is given byi 



This mapping is such that the point at infinity in the physical 




The mapping function isi 



(56) 





plane is mapped 




flow is then in the upper half of 



the £ plane. Figure 4(a) shows the resulting plane. 
The solution for the complex velocity isi 



(5 7) 



w(£) _ + o=" 

^ V£(£- b)' 



where , 
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and 




With equation (57) is substituted into equations (21) 



and (22) over the cavity and wetted surface of the foil, the 
vortex and source distributions can be obtained* Figures 5 
and 6 show the vortex and source distributions for cavity equal 
to one-half chord length. These figures show the leading edge 
and cavity termination singularities* One of the difficulties 
for the numerical procedure is to reasonably approximate these 
four singularities* 

Other quantities of interest are the ratio J2L , the 
lift coefficient and the cavity area* These are given ast 






-/ + 3 sinS) cos S sin $ + 

6 (l+sin£)(i+ 3sih $~Zsin 2 £- 6sin 3 S)^ 



+ -J- coian 



and 
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The equation for cavity area given in £ l] is in error 
and was corrected to equation (60) in [2J « Plots of equations 
(58) through ( 60 ) are shown in figures 7» 8 and 9* 

To obtain a solution to the super cavitating problem, 
Geurst performs a similar mapping (figure 4(b) )defined byi 





The solution for the complex velocity is then. 




where 




jf " cos *<$ 



(F 



JHEKL \lar>S - CJ>sC%-4z)} 

3 s 



and 



B _ fs in S' -tcu n & Cos ~~ ) ~f~ SI *1 ~jt)J 
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The vortex and source distributions are obtained using 
the same procedure as with the partial cavitation case.. Figures 
10 and 11 depict these distributions. In the super cavitation 
case only a leading edge singularity appears in the vortex 
distribution. Again the source distribution shows the leading 
edge and cavity termination singularities. 

The quantities ~r » lift coefficient and cavity area 

are « 

-f = i 



(6+) 



Cl 



/ 

; % sin SC 1 + smS) 



and 



jifr cc 



e<p s & 

3 H sinSO-s/nS)* 



Cbb) 
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Figures 12, 13 and l4 depict equations ( 63 ) through (65)* 



As with partial cavitation when the cavity approaches chord 
length the lift coefficient and cavity area become singular. 

In chapter VI the computed results are compared to the 



analytic solution* 



CHAPTER VI 
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RESULTS, CONCLUSIONS AND RECOMMENDATIONS 

The primary result of this investigation is a numerical 
method for two-dimensional cavitating flows which gives computed 
values close to the analytic solution. The computed results 
have been plotted in figures 3 and 5 through 16. In figures 5 
through 14, the analytic results of Geurst are also shown. 

Figures 5 and 6 show the computed vortex and source 
distributions for a cavity of one-half chord length. Each com- 
puted point is plotted at the control point location, except 
for the last vortex element which is plotted at x^ M • It is to 
be remembered that each computed point represents a constant 
value from Xg £ to %j iH . These two figures show that the com- 
puted distributions fit the analytic solution reasonably well. 
However, the largest magnitude elements are not shown on the 
figures (see the sample computer output in Appendix B for their 
values). These elements do not appear to approximate the 
proper singular behavior. 

Near the singularities the distributions of sources and 
vorticies can be approximated by l/x n . For the leading edge 
n has the value one-quarter and for cavity termination n 
has the value one-half Cl]. The computed values of n are 
given in Table 4 for three values of M and N. From Table 4 it 
is apparent that the behavior of the two closest elements does 
not match closely the analytic solution's singular behavior. 



Figures 7» 8 and 9 compare the computed values of , 

(■‘■/in'* and A/ifrot. * for ten cavity lengths to the analytic 
solution* The plotted results show the points closest to the 
analytic solution from a number of computations for various 
values of M and N. In general the computed results compare 
quite favorably with the analytic solution. The exception to 
this is for cavity lengths approaching chord length. For these 
cavities the method appears to be very slowly convergent. For 
these cavities a large number of elements are required for the 
cavity-foil region leading to a large number of equations. 

Table 3 contains all the computed results to date. This 
table gives results for variations of M and N for specific 
cavity lengths. For example, taking X 9 0.5» the results show 
that when the vortex and source element size is smaller in the 
combined cavity-foil region than in the fully wetted region, 
the computed values show an improved level of convergence to 
the analytic solution. This result is, also shown in figures 
5, 6, and 15 * 

Figure 15 shows the results of a convergence test for a 
one-half chord length cavity • this figure also shows the 
tendancy for better convergence for a smaller interior element 
size. It is noted that for figure 15 the vortex element size in 
the region of fully wetted flow are the same for both ratios 
of M and N. This effect is also shown in figure 5 for the vortex 
distribution with the same ratios of M to N and the source 
distribution in figure 6. 



This conclusion points out one of the problems with 
the numerical method used herein. The method reported on require 
that the vortex and source elements in the combined cavity- 
foil region be of equal size. An independent variation of the 
source and vortex element sizes cannot be performed. Thus it is 
not possible to determine whether decreasing the size of the 
vortex elements, the source elements .or both, results in the 
improved convergence. It is recommended that an investigation 
of effects of independent variations of the vortex and source 7 
elements be done. This can be achieved through a moderate re- 
vision of the computer program. 

Also, since the singularities do not appear to be well 
matched, a further investigation of the location of control 
points should be considered. It seems desirable to have an 
analytic basis for the control point locations. 

From the general form of the computed vortex and source 
distributions, this investigator is convinced that asstiming 
piecewise linear distributions will produce an immediate improve 
ment in the computational results. It is, therefore, recommended 
that this numerical method be modified to incorporate piece- 
wise linear vortex and source distributions. 

The results for the super cavitation case are very 
similar to the partial cavitation case. The convergence to the 
analytic solution becomes worse as the cavity length approach 
chord length. The same conclusions and recommendations made 
for the partial cavity apply to the super cavity. 
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Another purpose of this investigation has been to 
develop a method of allowing arbitrary values of angle of attack 
and camber as inputs and then determine the cavity length. An’ 
itterative procedure to accomplish this has been developed. 

The method and the computer program are in Appendix C with the 
computed results. The method and program are a preliminary effort 
provided to prove only that such a method is possible. 

In summary these are five recommendations for future 
investigation. 

1 ) A program be written that allows the vortex 
and source element sizes to be determined 
independently. 

2) An analytic and further numerical investigation 
of control point locations. 

3) Piecewise linear distributions of sources 
and vorticies should be investigated. 

4) That the method herein not be applied to the 
three-dimensional case without further 
investigation. 
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5) The program should be modified to include 
cambered profiles. In this case the ratio 
of <*/CT cannot be computed as an unknown. 
Either the angle of attack or cavitation 
number must be known along with the camber 
line. It is not clear that the control 
point locations will be the same as with 
the flat plate. 

6) Further development of the method of 
computing cavity length from arbitrary 
values of <x and c r • 

The conclusions drawn from this investigation are« 

1) the method gives computed values close to 
the analytic solution. 

2) the distribution of sources and vorticies 

is probably the best test of computed results, 

3) the method shows poor convergence for 
cavity lengths near chord length. 
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4) an iterative procedure to determine cavity 
length from an input of angle of attack 
and cavitation number is possible* 
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figure 6 

Source Distributions J?=0.5 
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figure 8 

Lift Coefficient - Partial Cavitation 
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Table 1 

Values of M and N for figures 5 through 9 



Cavity length 
Si 


No. of Vortex 
elements 
M 


No. of Source 
elements 
N 


0.1 


30 


15 


0.2 


30 


15 


0.3 


30 


15 


0.4 


30 


15 


0.5 


54 


34 


0.6 


25 


20 


0.7 


25 


20 


0.75 


30 


25 


0.8 


30 


25 


0.9 


30 


25 


0.95 


30 


25 



Geurst 

Computed M=IO N=20 




Vortex Distribution , X = l-5 
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figure 1 1 Source D is trbutior). j ?=/,5 
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figure 13 



Lift Coefficient -Super Cavitation 
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Cavity Area-Super Cavitation 
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Table 2 

Value 8 of M and N for figures 10 through 14 



Cavity length 

9 


No. of Vortex 
elements 
M 


No. of Source 
elements 
N 


1.05 


20 


25 


1.10 


20 


25 


1.50 


15 


30 


1.60 


15 


30 


1.80 


15 


30 


2.00 


15 


30 


2.50 


15 


30 


3.00 


15 


30 


3.50 


15 


30 


4.00 


15 


30 


5.00 


15 


30 
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Numerical Convergence, j 2=05 
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figure 16 
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Table 3 

Computed Results 
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cavity 


No. of 


No. of 


(X 


Cc 


Area. 


elapsed 

time 


length 


Vortex 


source 


cr 


' Sl'TUX 


oi 7iC< 




elements 

M 


elements 

N 




* 

.1 


30 


15 


.0380 


1.039 


4. 66x10“ 3 
3. 45x10“ 


40 


.1 


45 


5 


.0485 


1.002 


57 


.2 


30 


15 


.0569 


1.066 


.0129 


41 


• 2 * 


4o 


10 


.0627 


1.035 


.0119 


61 


.3 


30 


15 


.0718 


1.0913 


.0237 


42 


.4* 


30 


15 


.0844 


1.12 


.0368 


40 


.5 


10 


5 


.107 


1.09 


.0412 


2 


.5 


20 


10 


.0984 


1.14 


.0495 


15 


.5 


30 


15 


.0953 


1.15 


.0520 


42 


.5 


40 


20 


.0936 


1.16 


.0533 


94 


.5 


62 


31 


.0916 


1.18 


.0549 


324 


.5 


30 


25 


.0796 


1.33 


.0590 


72 


•5* 


27 


17 


.0913 


1.19 


.0536 


38 


.5 


54 


34 


.0891 


1.20 


.0558 


283 


.5 


25 


20 


.0819 


1.30 


.0576 


42 


.6 


30 


25 


.0882 


1.39 


.0785 


70 


.6* 


30 


15 


.1047 


1.20 


.0695 


43 


.6* 


25 


20 


.0909 


1.36 


.0765 


4o 


.7 


30 


25 


.0952 


1.48 


.101 


70 


.7* 


30 


15 


.112 


1.26 


.0898 


42 


.7* 


25 


20 


.0982 


1.44 


.0984 


40 


• 7 $ 

• 8* 


38 


\\ 


m 


1.537 

1.61 


.114 

.129 


68 

73 




25 


20 


.103 


1.56 


.125 


43 


.9 


30 


25 


.0981 


1.88 


.167 


73 


.9 * 


25 


20 


.102 


1.81 


.161 


41 


.95 


30 


25 


.0906 


2.22 


.199 


71 


.95 


25 


20 


.0947 


2.12 


.191 


40 



* 



Points plotted in the figures. 



Table 3 continued 



cavity 

length 


No. of 
Vortex 
elements 
M 


No. of 
source 
elements 
N 






fi ^ rca . 

an* 


elapsed 

time 

( 7 To S « c ) 


cr 




1.0$* 


20 


25 


.136 


1.55 


.192 


41 


1.1 


20 


25 


.185 


1.082 


.177 


44 


1.5 


10 


20 


.384 


0.506 


.230 


15 


1.5 


5 


10 


.413 


0.445 


.217 


3 


1.5* 


15 


30 


.379 


0.515 


.233 


39 


1.5 


20 


40 


.372 


0.527 


.236 


89 


1.5 


17 


27 


.379 


0.514 


.234 


38 


1.5* 


20 


25 


.405 


0.461 


.226 


40 


1.6* 


15 


30 


.414 


0.475 


.254 


42 


1.8* 


15 


30 


.477 


0.422 


.300 


41 


2.0* 


15 


30 


.532 


0.388 


.348 


43 


2.5* 


15 


30 


.653 


0.338 


.482 


41 


3.0* 


15 


30 


.757 


0.311 


.620 


40 


3.5* 


15 


30 


.851 


0.293 


.790 


41 


4.0* 


15 


30 


.937 


0.281 


.960 


40 


5.0 


15 


30 


1.093 


0.263 


1.330 


42 



65 



Table 4 

Singular Behavior 



M 


N 


n* 

leading 

edge 

Vortex 


n* 

cavity 

termination 

Vortex 


n* 

leading 

edge 

source 


n* 

cavity 

termination 

source 


20 


10 


•093 


3U1 


-.794 


36 


27 


17 


.099 


27.5 


-.716 


35.5 


60 


30 


-.106 


-60.5 


-.604 


57.7 



* Based on only the two elements closest to 
the singularity (e.g. for the leading edge 
vortex distribution the elements are the 
first to vortex densities) and j? * 0*5* 
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APPENDIX A 



COMPUTER PROGRAM 



Appendix A includes the flow charts for the computer 
program and the computer program listing. The flow charts 
are provided to give a broad overview of the program flow. 

The numbers between blocks in the flow charts indicate 
approximate statement numbers in the program listing. 

A flow chart for RMINV is not provided since it is a standard 
matrix inversion routine. 
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DATA ND,NDIK/100, 101/ 

DATA IN , OUT, IPCH/5 , 6 , 7/ 
PI=2.*ARCOS (0.0) 

READ (IN , 1 000) H,N,IWRT,CAVL 
IF (H .EQ. 999) STOP 
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25 CONTINUE FMNV0320 

RHNV0330 

INTERCHANGE RONS. RMNV0340 

I = I WORK (K, 1) RMNV0350 

IF (I .LE. K) GO TO 35 RNNV0360 
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DO 55 I = 1, NORDER RMNV0680 

IP (I . EQ. K) GO TO 55 RMNV0690 

HOLD = - A (I, K) / BIGA RHNV0700 

A (I, K) = HOLD RHNV0710 

RHNV0720 
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130 CONTINUE RBNV1040 

140 CONTINUE RMNV1050 

RfINV 1060 

IERB = 0 RMNV1070 

RETURN RMNV1080 



RMNV1090 

1001 FORMAT (23H RMINV: ARGUMENT ERROR, 2111) RMNV1109 

1002 FORMAT (28H RMINV: MATRIX IS SINGULAR.) FMNV1110 

RMNV1120 

END RMNV1130 
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APPENDIX B 



SAMPLE COMPUTER OUTPUT 



Appendix B contains the sample computer output. This 
output contains the results of a convergence test for 
cavity lengths equal to £ chord length and 1.5 times chord 
length. Note that for a partial cavity source control 
points beyond cavity termination are listed in the output. 
These control points are superfluous and not used in the 
program. The same is true of control points for the 
vortex elements beyond chord length in super cavitation. 

The listing of source and vortex densities has the 
same order as the listing of control points. 
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NO. CF V ORTICI ES= 1C 
NO. OF SO URCE S= 5 

CAVITY LENGTH = 0.50000 

*** DETERM IN ANT=-.1 317091 E+09*** 

LIFT COEF TC IE NT/ (2*PI*A LFA ) = 
AREA OF CAVITY/ (2*PI*ALPA) = 
TOTAL SOORCE STRENGTH/SIGMA = 
ANGLE OF ATTACK/SIGMA 
CAVITY LENGTH 

ELAPSED TI ME = 2 1/100 SEC. 



* 



0. 10 926 49E + 0 1 
C. 41 17717 E- 01 
C . 2272427 E-06 
0. 1067560E+00 
0. 5000000E+00 
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NO. OF VORTICIES 3 20 
NO. CF SOU EC ES= 1C 
CAVITY LENGTH= 0.50000 



XL(I) 

0.0 

0. 31 2 5 0E-01 
0. 6 2500E-0 1 
0.1 2500E+00 
0. 18750E + 00 
C. 2500 0E +00 
0. 3 1 250E+00 
0.37500E+00 
0. 4375 0E + 00 
0. 46875E+00 
0.50000E+00 
0. 53 1 25E+00 
0.56 250E+00 
0 . 6250 CE + 00 
0. 687 50E+ 00 
0.75000E+00 
0. 81250E + 00 
0. 8750 0E + 00 
0.93750E+00 
0. 96875E+00 

*** DETERMINANT 3 



X (I) 

0. 28125E-01 
0.59375E-01 
0. 1 1875E+00 
0.18125E+00 
0. 24375E+00 
0.30625E+0G 
0. 36875E+00 
0. 43125E+00 
0.46562 E+00 
0.49687E+00 
0.5281 2E+00 
0.55937E+00 
0.61 875E+00 
0.68125E+00 
0.74375E+00 
0.80625E+00 
0.86 875 E+00 
0.9 3 125E+00 
0 . 96562 E+00 
0. 99687 E+00 

. 3532951E+ 17** 



XU(I) 

0. 31250E-01 
0.62500E-01 
C. 12500E+00 
0. 18750E+00 
C.25000E+00 
0.31 250E+00 
C. 37500E+00 
0. 43750E+ 00 
0 .46 875E + 00 
C. 50000E+00 
0. 53 1 25E+00 
C.56250E+00 
C.62500E+00 
0.68750E+00 
C.75000E+00 
0.81 250E+00 
0.87500E+00 
0.93750E+00 
0.96875E+00 
C . 1 0000E+0 1 



VORTEX DENSITIES/SIGMA: 



-0.93763E+00 
- 0. 6856 CE+00 
-0. 64077E+0C 
-0.45946E-02 
-0. 63426E-01 



-0.87444E+00 
-0.66212E+00 
-0.55656E+00 
-0.59571 E-01 
-0. 51625E-01 



-0.79342 E+00 
-C.63783E+00 
0. 65990E+00 
-0. 71 848E-0 1 
-C. 37786 E-01 



SOURCE DENSITIES/SIGMA: 

0. 805 1 5E+00 0.53662E+00 

0. 10949E+00 0.76209E-02 

-0.56007E + 00 - C. 1 3 084 E +01 



C. 33202 E+00 
— C. 89 1 98E- 0 1 



LIFT COEFTCIENT/ (2*PI*ALFA) = 0. 1 1 37383 E+01 
AREA CF CAVITY/ (2*PI*ALFA) = C . 4946671 E-0 1 



TOTAL SOURCE STRENGTH/SIGMA =0. 1750886E-06 



ANGLE OF ATTACK/STGM A =C. 9841657E-01 



CAVITY LENGTH 



=C . 5000000E+00 



XS (I) 

0. 28125E-01 
0. 46875E-0 1 
0. 93750E-01 
0. 15625E+00 
0. 21875E+00 
0. 28125E+00 
0 • 343 75E+00 
0. 40625E+00 
0. 453 1 3E+00 
0.47187E+00 
0. 51563E+00 
0 . 54688E+00 
0. 59375E+00 
0. 65625E+00 
0.71875 E+00 
0. 781 25E+00 
0. 84375E+00 
0.90625E+00 
0. 95313E+00 
0 .98438 E+00 



0. 72474E+00 
0. 6526 1E+00 
0. 11008E+00 
0.7 09 11 E-01 
0. 38749E-01 



0. 20736E+00 
0. 30393E+00 



ELAPSED TIME 3 



15 1/100 SEC 
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NC. CF VORTICIES= 10 
NO. OF SO URCES= 20 
CAVITY LENGTH= 1.50000 



XL {I) 


X (I) 


XO(I) 


XS(I) 


0.0 


0.56250E-01 


0.62500E-01 


0.56250E-01 


0. 6 250 OE-O 1 


0. 1 1875E+00 


0. 12500E+00 


0. 93750 S-01 


C. 1 2503E+00 


0.23750E +00 


0.25000E+00 


0. 18750E+00 


0. 25000E+00 


0.36 250E + 00 


C. 37500E+00 


0. 31250E+00 


0. 37500E+00 


0.48750E+00 


0.50000E+00 


0. 43750E+00 


0. 50000E+00 


0. 61 250E+00 


C.62500E+00 


0 . 56250E+00 


0.62500E+00 


0. 73750E+0C 


C.75000E+00 


0. 68750E+00 


C.7500CE+00 


0.86250E+00 


0 . 87500E+00 


0. 81250E+00 


0. 87500E + 00 


0. 93125E+00 


C.93750E+00 


0. 90625E+00 


0. 9375 0E+00 


0.99375E+00 


0. 10000E+0 1 


0. 96875E+00 


0. 1 OOOOE+OI 


0.10281 E+01 


0.10 3 13E+0 1 


0. 10 156E+0 1 


0. 10313E + 0 1 


0. 10 594E+0 1 


C. 106251+01 


0. 1 0469E + 01 


0.10625E+01 


0.111 87E+01 


0.11 250E+0 1 


0. 10938E+01 


0. 1 1250E + 0 1 


0. 1 1 81 2 E+01 


0.11875 E+01 


0.1 1563E+01 


0. 11875E+0 1 


0. 12437E+01 


0. 12500E+0 1 


0.121 88E+01 


C .1 2500E+01 


0.13 062 E+01 


0.13125E+01 


0. 1281 3E+01 


0. 13125E + 01 


0. 13687E+01 


C. 13750E+01 


0. 13438E+01 


0. 13750E+01 


0.14312S+01 


0. 14375E+ 0 1 


0. 14063E + 01 


0. 14375E+01 


0. 14656E+01 


0. 14688E+0 1 


0. 14531E + 01 


0. 14688E + 01 


0. 14969E+0 1 


C. 15000E+01 


0. 1471 9E+01 


*** DETERMINANTS 


. 3851246E+13*** 







VORTEX DSNSITIE; 


S/SIGMA: 






-0. 11499F+01 


-0. 10381 E+01 


-C. 87268F+00 


-0.72748E+00 


-0.62812E+00 


-0. 54372E+00 


- C. 460 15E+00 • 


-0.37014E+00 


-0. 26088 F. + CO 


- C, 25544 E+00 






SOURCE DENSITIES/SIGMA: 






0.1 3267F+0 1 


0.96 160E+00 


0.7 3057E+00 


0.6011 2E+00 


0. 50964E+00 


C.43215E+00 


C .36 176E+00 


0 . 285 1 OE+OO 


0. 2508 4E + 00 


0. 23408E+00 


0. 19550E-04 


-0. 33281E-01 


-C.164S1E+00 


-0. 37803 E+00 


-0.49625E+00 


-0.75565E+00 


-0. 9853 2E + 00 


-0. 15239E+01 


-C.29826E+01 


-0.56046E+01 



LIFT COEFICIENT/ (2*PI*ALFA) =0.5062250E+00 
AREA OF CAVITY/ (2*PI*ALFA) =C . 2296074 E+00 
TOTAL SCUPCE STRENGTH/SIGMA =- . 1788 1 39 E-06 
ANGLE OF ATTACK/SIGMA =0 . 3843980E+00 
CAVITY LENGTH =0.15030003+01 



ELAPSED TIME= 



15 1/100 SEC 
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NC. OF V ORTICI ES= 15 
NO. OF SOURCES^ 30 
CAVITY LE NGTH= 1.50000 

*** DETERMINANT=0.7520241E+19*** 

LIFT COEF ICIENT/ (2*PI*ALFA) = 
AREA OF CAVITY/ (2*PI*ALFA) = 
TOTAL SOURCE STRENGTH /SIGH A = 
ANGLE OF ATTACK/SIGMA 
CAVITY LFNGTH 

ELAPSED TIME= 39 1/100 SEC. 



0. 5145274E+00 
C.2330343E+00 
C. 1430511E-05 
0. 3788325E+00 
0. 1500000E+01 



APPENDIX C 



COMPUTER PROGRAM FOR COMPUTATION 

OF 

CAVITY LENGTH FROM ARBITRARY VALUES 

OF 



ANGLE OF ATTACK AND CAVITATION NUMBER 



95 

The iterative method for computing cavity length from 



arbitrary values of angle of attack and cavitation number is 
based on assuming an initial cavity length ( J? * 0.5 ) then 
computing the next cavity length ast 



In equation (1 C) the value of (■§=); is the value computed 
from Jt L . This works well for moderate length partial cavities. 

For super cavities the method used is to set upper 
and lower boundaries on the cavity length. Then the next 
assumed value of cavity length is* 



determined by comparing to input. For ($L) input 

greater than [~r) L the actual cavity length is greater than Jii 
and so long as / L is greater than/^the lower boundary for the 

cavity length is then • For less than (%.) c 

Xi becomes the new upper boundary. Thus the boundaries on 
cavity length will converge to a solution. 

The following program listing uses this procedure. This 
program also uses the subroutines CPGEN, MATRIX and RMINV. 



(1C) 





The lower and upper boundaries on cavity length are 
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The iterative procedure for the super cavitating case 
was developed to overcome apparent nonconvergence of the method 
used in partial cavitation# However, the resulting nonconver- 
gence was based on theoretical cavity lengths only slightly 
greater than chord length (e.g# Jl ^ 1.04). In this region 
of cavity length <*/c r increases very rapidly with cavity 
length# Convergence of any method for these cavity lengths 
will be slow. Therefore, the method used herein is probably 
an artificial result and not a required procedure. 
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Table C 1 

Results of Cavity Length Computations 



input 

c* 

c r 


computed 

JL 

cr 


computed 

cavity 

length 


No. of 
iterations 


Geurst 

cavity 

length 


Percent 

difference 


.025 


.0249 


• 0444 


4 


.045 


-1-3 


.05 


.0498 


.145 


3 


.155 


-6.5 


.075 


.0745 


.311 


3 


.375 


-17.1 


.10 


.998 


.730 


3 


1.04 


-29.8 


.15 


.1498 


1.06 


6 


1.099 


-3.5 


.20 


.199 


1.12 


9 


1.16 




.25 


.249 


1.19 


6 


1.25 


- 4.8 


.333 


.335 


1.34 


6 


1.45 


-7.6 


.50 


.498 


1.875 


4 


2.0 


-6.25 


1.0 


.9998 


4.39 


8 


5.0 


-12.0 
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2000 FORMAT (» 1***ERR0R IC0UNT= 10**** , 3E20 . 5) 
2010 FORMAT('0' ,15X,3E15.5) 

2020 FORMAT ( ' 1 ' ) 

2030 F0RMAT(‘ »,E15.5) 

END 
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